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The temperature dependence of the frequencies of a Bose-Einstein condensate obtained 
in experiment has not been fully understood theoretically. In this paper we present a 
simplified version of a two-gas model. A numerically-found ground state of the system is 
used for the small-oscillations analysis. In the case of spherical symmetry a full spectrum 
of frequencies is found for low orbital quantum numbers. Avoided crossings that appear 
in the spectrum might be the reason for experimentally observed frequency shifts. 

I. INTRODUCTION 

One of the first experiments done after achieving Bose-Einstein condensates was measuring its eigen- 
frequencies . These experiments were performed at very low temperatures, where no thermal atoms 
were detectable. The results agreed very well with the predictions based on the mean-field theory |^ ||. 
This success was not repeated with later experiments, which measured the frequencies of the conden- 
sate in full temperature range Jt],||]. To explain the experimental results several models were proposed 
P~H • To the best of our knowledge none of these models managed to reproduce experimentally observed 
frequency shifts. 

In this paper we investigate a simplified version of a two-gas model. The thermal cloud influence on 
frequency spectrum is taken into account, as well as dynamical interaction between the gases. Thermal 
atoms are treated in hydrodynamic limit but no mechanism for exchanging the particles between the 
phases is provided. A ground state of such a system is found numerically, and it shows, that there 
are qualitative differences with a ground state considered in other models. For small oscillations linear 
response theory is applicable - it leads to an eigenproblem. In the case of spherically symmetric trap the 
resulting equations are solved numerically. 

The paper is organized as follows. In Section 2 the model is furthermore specified. Numerical strategies 
are described in Section 3, while Section 4 is devoted to the results. 



II. TWO-GAS MODEL OF BEC AT FINITE TEMPERATURE 

We develop a simple model of Bose-Einstein condensates at finite temperatures. In this picture the 
trapped Bose gas at temperature below the Bose-Einstein condensation temperature is considered as 
a system consisting of two distinct components. One of these components, the condensed atoms, is 
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described by the generalized Gross-Pitacvskii equation. The second component, the cloud of thermal 
atoms, behaves like a classical ideal gas. We assume that the relaxation processes within the thermal gas 
are extremely fast. Hence, the thermal gas is all the time in the thermal equilibrium and its temperature 
can be determined experimentally by means of its free expansion. The model works in the so called 
hydrodynamic regime as opposed to the approach presented in [ fl3| investigating collisionless regime of 
parameters. Similar model was considered in [Q. 

The number of atoms in each component is preserved, i.e., no mechanism allowing the exchange of 
atoms between two factions is included in this version of the model. The fluctuations of the number of 
atoms in the condensate have been calculated to scale with the number of atoms as y/Afo jl5| ■ Therefore 
their influence is negligible (less than 1% for temperatures below 0.9 T c in our simulations) since we are 
dealing with J\f—50 000 atoms. The relative number of condensed and thermal atoms is taken from the 
thermodynamic expression for an ideal gas in the trap. As we know (HQ, residual interactions and 
finite size effects modify the temperature depletion curve 1 — (T/T c ) 3 in a qualitatively insignificant way. 

The model is time-dependent. It is defined by the following set of equations: 
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where ^(r, t), pth(?, t), and v^(r, t) are the condensate wave function, the density and velocity fields of 
the thermal cloud, respectively. 

Contrary to the models which consider the thermal part of the Bose gas as a static cloud [ ^2Ul^Jl9] ] , 
in our approach both components dynamically affect each other. The influence of the thermal part on 
the condensed one is described by a kind of Gross-Pitaevskii equation with an extra term (the last term 
on the right hand side of the first equation of (|l|)). This term can be considered as a time-dependent 
generalization of the mean-field interaction with the non-condensed density derived in Hartree-Fock- 
Bogoliubov-Popov approximation pofl . As has been shown in jL2|, going beyond Popov approximation 
by taking anomalous averages into account does not reproduce experimental results. Hence, in this paper 
we concentrate on different mechanisms that can lead to frequency shifts. The second and the third 
equations of ([j]) describe the thermal cloud as a classical gas confined in the effective potential formed 
by the external potential and the condensate density (the last term in the third equation of ([l])). 

In the case of no coupling between two components we have independent oscillations of pure condensate 
(in accordance with the Gross-Pitaevskii equation) as well as the thermal cloud (in a hydrodynamic limit). 
Assuming no macroscopic motion of the thermal component, the third equation leads to the Maxwell- 
Boltzmann distribution for a gas confined in the external potential. 

Before going into further analysis of the equation ([!]), the trap potential should be specified as an 
anisotropic harmonic oscillator: 
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V ext = — (oj z x x z + u z y z +uj z z z J } 

Then the external potential defines both the length scale (the size of the ground state wavefunction) 
and the time scale r = u>q , where luq is the geometrical average of trapping frequencies 



ujq = ^LO x uj y uj z . Then the equations ([!]) can be expressed in the following form: 
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Pth{r, t) 
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nin m (r,i) +E^ + 9 v th( r >t) + 2 7 Pc(r,t) 



Now LOi = Lu Xi /ujQ, denotes the relative temperature: f2 = 2kT/fiu)o, and 7 measures the strength of 
interactions: 7 = 16na (a is dimensionless scattering length). 

Stationary solutions of (|^) are found by using the imaginary-time propagation technique. First, the 
second and third equations of (^) are rewritten in terms of a nonlinear Schrodinger equation, which is 
allowed provided the thermal velocity field is irrotational (this assumption was actually utilized in writing 
the third equation of (|])). An extra term appears on the right hand side of the second equation of (|J); 
it is up to the sign the so called "quantum pressure" term, known from the hydrodynamic formulation 
of quantum mechanics fell. 
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To find the stationary solutions of (||) one has to make transition to the imaginary time: t — > i t. 
Assuming, the ground state wave function of a two-gas system described by (|[) is real, we obtain the 
set of equations for real functions $ c (f, t) and $t?i(r, t) which in the limit t — > 00 converge to the ground 
state wave functions: V'cC 1 ") an d V'tftC 1 ')- 
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FIG. 1. Densities of the ground state in the case of spherical trapping potential. Solutions of (^) are plotted 

with continuous lines, dotted lines denote the ground state of the system in the case of no interaction between 

phases. On the left figure there is 5 000 atoms, on the right 5 000 000. On both the temperature is 0.67 T c . 

Examples of the ground state densities are plotted in figure @. If there was no interaction between the 
phases, we would get a parabola for the condensate and a Gaussian for the thermal cloud. Taking this 
interaction into account, however, results in slight change of the condensates density, and in great change 
of the thermal cloud, which is no longer described by the Gaussian, but by a function with a maximum 
at the edge of the condensate. 

The imaginary-time propagation technique was proved to work in the case of linear Schrodinger equa- 
tion. We have checked that it also works for the highly nonlinear problem just defined. The set of partial 
differential equations (||) constitutes an initial-value problem of the parabolic type. We have solved it by 
using various numerical methods, depending on the dimensionality of the system (for a one-dimensional 
Bose gas the split-operator technique was used whereas in a three-dimensional case we applied the 
Runge-Kutta algorithm) . 

To find the spectrum of frequencies in the two-gas model, we use the method of small oscillations 
around the numerically found ground state of the system. The assumption of small oscillations can be 
written: 



p c {r,t) = p°(r)+. 



<5p c (r) 



Pth(T,t) = p Q th {T) + e-^Sp th (v) 
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where the common frequency for both phases is already introduced. Of course, p° c (r) = (V£(r)) and 
Pt h (r) = (V'tk( r )) 2 denote densities of the ground state of the condensate and the thermal cloud, re- 
spectively. Substitution these equations into the last two of (|3|), linearization in Sp c (r) and 8pth{^), and 
elimination of velocity field yields the following: 
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It is worth noticing that the equation (Q) does not allow for an arbitrary phase shift between oscillations 
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of the condensate and the thermal cloud. This can be easily seen by writing the equations (Q) in real 
(not complex) form and following the linearization procedure again. The resulting equation will have a 
time-dependent factor in it; setting this factor constant will lead to allowed phase shift: <f> = or (j) = iz. 
This argumentation is valid locally, which means that for high frequency modes in some spatial region 
both gases might oscillate in phase, while in others out of phase. 

As far as condensate is concerned, the small oscillation analysis leads to the Bogoliubov equations, 
modified by the presence of the thermal cloud. First, the second of the (^) together with the following 
ansatz: 



*(r,t) 



= e -iM*/?» 



(vO\W>°(r) + U (r)e- iwt + i;*(r)e iwt ) (8) 

are inserted into the first of the equations (||). Linearization procedure and comparing the coefficient 
standing next to e ±luJt results in: 



-Cu(r) + %p°(r)v(r) + 7 V /p°<5pth(r) = utt(r) 



(9) 



where C := -V 2 + ^ + 7 (p°(r) + ^(r)) - p. 

It is easily seen, though, that the role of 5p(r) in equation (^) plays \J p®(r) (u(r)+v(r)). The Bogoliubov 
equations can be transformed so to introduce 5p c := u + v: 



£ - Ip°(r)) { (L + 2pP( r )) Sp c (r) - 2 7a /^W Sp th (v)} = u? 5 Pc (v) 
The equations (ffl) and (HQ) constitute an eigenvalue problem: 
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nv-(p°„v4-) 



where C := -V 2 + £. ^ + 7 (p°(r) + - p. 

Up to now the considerations were fully anisotropic. For further analysis however, only the simplest 



case, the spherical symmetry is assumed. Then the external trap is defined: V ex t = 



In systems 



with spherical symmetry the angular dependence takes form of the spherical harmonics; it is then possible 
to seek the global solution in the form: 



Sp c (r) = F(r)Y lm (6,<j>) 
5 Pth {r) = Q{r)Y lm {0,<t>) 



Then the equation (p. l|) transforms to: 
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III. NUMERICAL APPROACH 



The problem is first discretized on a lattice. Then the equations ( |l3| ) are solved directly by diagonal- 
ization of the resulting matrix. The full spectrum of eigenfrequencies and eigenmodes is found. 

Boundary conditions are set both at the boundary of the lattice and at the origin. For large distances 
we demand that the functions T{r), •F'(r) and Q{r) vanish. In the origin the boundary conditions depend 
on the value of orbital quantum number I. 

To obtain these we used the analytical expression for the condensates wavefunction valid in the Thomas- 
Fermi limit and in T — regime || : 



jF(r) = V4n + 21 + 3 yffc r l p( l+1 ^ (1 - 2r 2 ) (14) 



where p^ +1 / 2 ' ) are Jacobi polynomials fl22fl . For I = the boundary conditions are: T 1 (r = 0) = and 
T'"(r = 0) = 0; for I = 1: T{r = 0) = and T" (r = 0) = 0. Of course, they are identical as for the 
radial equation of the ordinary linear Schrodinger equation. 

For the thermal cloud, the equation (0) in the lack of interaction between the gases can be solved 
analytically. The ground state is just a Gaussian: 

&(r) = 7S^«p(-£) (15) 



(27rfl)3/ 2 ^ V 2^ 
and the second of the equations ( |l3| ) takes the form: 

ng"(r)+g\r)(—+r)+g(r)(u 2 + 3- l ^fin)=0 (](>) 



r 



Two independent solutions are: 



(17) 



where 1F1 denotes Kummer's hypergeometric function p2fl . The solution Q^ 2 \r) is singular at r = 0. For 
the solution Q^{r) to vanish at r — > oo the frequency must be: 



uj = V2n + I (18) 

From the explicit form of the solution the boundary conditions can be deduced: G'{r = 0) = for / = 
and G(r = 0) = for I > 0, again as in quantum mechanics. 

IV. RESULTS 

All of the calculations were performed for the rubidium condensate with parameters: m = 1.44- 10 _22 g, 
a = 5.821nm, uiq = 2tt ■ 200Hz and M = 50 000 atoms. This does not mean however, that the results 
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are specific to this set of parameters only. In our model there are only two dimensionless parameters 
describing the system: 7 and 17. The type of atoms comes into play only when 7 ~ a/y/m is concerned. 
For the sodium and rubidium condensates the difference in 7 is around 8%. 

Figures I and I are the main results of this paper. They show temperature dependence of the low- 
energy excitation spectrum for different values of orbital quantum number. For T = and T = T c 
frequencies are known analytically: dots represent the spectrum g: 

lj = \Jln 2 + 2nl + 3n + l (19) 

and squares the spectrum (fL8|). It can be seen that numerical results are in fair agreement with these 
spectra. The "horizontal" lines are the frequencies where thermal cloud oscillation is non-zero at all 
temperatures; for T — where there is no thermal cloud, those frequencies are not valid anymore. A 
complementary set of lines (inclined) originates from the T — condensate frequencies. 
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FIG. 2. Eigenfrequencies versus temperature for 1 = 0. 
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FIG. 3. Eigenfrequencies versus temperature for 1 = 1. 

For the temperatures between T = and T = T c both phases oscillate with similar amplitudes. 
Examples are plotted in figure 0. 
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FIG. 4. Shapes of eigenmodes. 



condensate 
thermal cloud 



10 15 20 25 30 
distance 



The existence of the avoided crossings in the spectrum in figures |^ and [| is clear. They occur when 
two eigenmodes as a function of a continuous parameter approach each other, and after a certain point 
they start to diverge. 

Avoided crossing is a well known phenomenon in physics. Its general properties are described by 
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the Landau-Zener theory. According to it, the character of eigenfunction "jumps" across the avoided 
crossing, as seen on the figure || If the parameter is time, two types of processes might be recognized: 
adiabatic (slow) and diabatic (quick). The time scale distinguishing between those types is connected to 
the minimal distance between the eigenenergies - via the uncertainty principle. 

It would be interesting to observe these avoided crossings experimentally Unfortunately, due to partic- 
ular realization of experiments, those processes cannot be observed, because temperature is a parameter 
given from outside, not changing during the experiment. However, our conjecture is that the observed 
rapid dependencies of the cigcnfrcqucncies with temperature occur because of the avoided crossings. 
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FIG. 5. Shapes of the eigenmodes around an avoided crossing. 

It is also clear that due to the interaction with the thermal cloud mode which in T = had fre- 
quency to = is awakening and no longer has it zero frequency. Therefore it has a potential of being 
experimentally observed. 
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V. SUMMARY AND ACKNOWLEDGMENTS 



We investigated a simplified version of the two-gas model. Within this model dynamic interactions 
between the phases: condensate and thermal cloud, are included. No mechanism for exchanging particles 
nor for damping is provided, though. 

Ground state of the system is found numerically using a propagation in imaginary time technique. 
Equations resulting from small oscillations theory constitute an eigenproblcm, which is numerically solved 
for spherically symmetric trap. 

The resulting spectrum of frequencies exhibits interesting features. Firstly, it is composed of lines 
originating from pure condensates oscillations in T = and from thermal cloud oscillations (T = T c ). 
In wide range of temperatures, however, both phases oscillate with approximately the same amplitudes. 
Secondly, avoided crossings are observed. We suggest that experimentally observed sudden shifts in 
frequencies might be a result of an avoided crossing. Finally, due to the interaction with the thermal 
cloud the awakening of mode n — is observed, and might be experimentally detected. To make contact 
with the existing experimental work we plan a similar study in the axially symmetric geometry. 

We would like to thank C.Clark for fruitful discussions. R.B. and K.R. acknowledge the support of the 
subsidy from the Foundation for Polish Science. M.B. and K.R. thank for the support of MCS Grant 
PAN/NIST-98-340. 
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